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Abstract.  We  derive  and  analytically  and  numerically  solve  statistical  moment  differential 
equations  for  immiscible  flow  in  porous  media  in  the  limit  of  zero  capillary  pressure,  with  application 
to  secondary  oil  recovery.  Closure  is  achieved  by  Taylor  expansion  of  the  fractional  flow  function  and 
a  perturbation  argument.  We  reduce  the  equations  by  exploiting  a  relationship  between  saturation 
and  velocity  correlations  that  is  unique  to  flow  in  one  dimension.  Mean  and  variance  of  (water) 
saturation  exhibit  a  bimodal  character;  two  shocks  replace  the  single  shock  front  evident  in  the 
classical  Buckley-Leverett  saturation  profile.  1 
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1.  Introduction.  Subsurface  geologic  properties  at  field  scales  are  uncertain, 
and  are  often  described  statistically  in  practice.  Flow  profiles  in  such  porous  media  are 
uncertain,  and  statistical  flow  outcomes  are  appropriate.  We  are  primarily  interested 
in  mean  behavior  and  a  measure  of  the  uncertainty  about  this  mean. 

A  “zeroth-order”  model  of  mean  flow  with  averages  of  geologic  properties  ignores 
correlations  between  flow  variables.  Monte  Carlo  simulations  of  many  realizations  of 
geologic  properties  to  estimate  moments  require  much  computation  time  and  careful 
sampling  techniques  [12,13,35].  Macrodispersion  theory  in  contaminant  transport 
captures  a  first-order  effect  of  fluctuations  via  covariance  functions  in  PDEs  for  the 
mean  concentration  [7, 14].  This  theory  has  a  long  history  in  subsurface  contaminant 
transport,  and  is  closely  related  to  eddy  diffusion  models  of  turbulence  [16,  p.  358  ff]. 

We  derive  second-order  PDEs  for  the  covariance  functions  and  the  mean  flow,  and 
solve  for  these  moments  simultaneously.  The  fundamental  problem  of  closure  of  the 
system  is  addressed  by  a  perturbation  argument.  The  resulting  moment  differential 
equations  (MDEs)  directly  approximate  the  local  mean  and  covariance  functions,  for 
general  boundary  conditions  and  general  stochastic  geology  [38]. 

1.1.  Applications.  A  statistical  description  of  subsurface  flow  is  of  particular 
interest  for  secondary  oil  recovery.  The  principal  difficulty  is  a  non-convex  nonlinear 
flux  function  in  an  advection  equation  that  leads  to  discontinuous  solutions.  Stan¬ 
dard  pressure-saturation  equations  for  2-D  horizontal  flow  of  two  immiscible  fluids  in 
porous  media,  in  the  limit  of  vanishing  capillary  pressure  ( pc  =  0),  are 

(j) v(x)  =  —  K(x)  Vh(x),  V  ■  v(x)  =  0,  (1.1) 

<9ts(x,  t)  +  V  ■  [/(s(x,  f))v(x)]  =  0.  (1.2) 

These  are  taken  to  be  valid  from  laboratory  (centimeters)  to  field  scales  of  reservoir 
depth  (10-100  meters)  and  length  (100-1000  meters).  Hydraulic  conductivity  K  may 
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be  an  anisotropic  tensor;  here,  for  simplicity,  it  will  be  an  isotropic  scalar  K.  Assume 
also  that  K  depends  weakly  on  (water)  saturation  s  [1,3].  Apply  (1.1)— (1.2)  to  the 
flow  of  oil  and  water,  for  arbitrary  fluid  mobilities.  Denote  the  total  velocity ,  a  scaled 
total  volumetric  flux  of  both  fluids,  by  v,  hydraulic  water  head  by  h,  and  porosity 
(assumed  constant)  by  <f>.  The  fractional  flow  function  f(s),  for  pc  =  0,  represents  the 
fraction  of  v  due  to  water.  It  is  typically  S-shaped,  and  we  use  a  form  arising  from 
quadratic  relative  permeabilities  for  illustration  (see  [1]).  However,  our  method  does 
not  depend  on  any  such  specific  choice  of  f(s). 

Capillary  pressure  regularizes  sharp  fronts  caused  by  the  nonlinear  advection 
term.  To  obtain  a  linear  approximation  to  this  effect,  add  ej^V2s(x,t)  with  en  >  0 
to  the  right  side  of  (1.2).  Letting  cd  0  defines  the  vanishing-viscosity  solution  [33], 
which  is  the  one  we  seek. 

As  is  standard  in  subsurface  applications,  let  Y  =  In  K  be  a  random  field  with 
prescribed  mean  and  covariance  functions;  e.g.,  it  is  often  claimed  that  Y  is  multi¬ 
variate  Gaussian,  based  on  empirical  observations  [10]  (our  method  does  not  depend 
on  this).  Through  (1.1)— (1.2),  v  and  s  are  thus  random  fields.  No  other  underlying 
sources  of  uncertainty  are  considered  in  this  study. 

Under  the  assumptions  stated  above,  with  steady  boundary  conditions,  and  ne¬ 
glecting  the  weak  dependence  of  K  on  s,  a  steady  v  can  be  determined  from  (1.1). 
We  evolve  s  from  the  stochastic  PDE  (1.2),  assuming  known  statistics  of  v.  Moments 
of  v  and  h  can  be  estimated  from  established  theory  ([36-38,41]  use  MDEs).  We 
combine  analytical  and  numerical  techniques  to  model  the  propagation  of  uncertainty 
from  an  underlying  random  field  Y (x),  through  v(x),  to  the  solution  s(x,t). 

1.2.  Previous  work.  Existing  work  on  MDEs  focuses  mostly  on  advection  equa¬ 
tions  with  linear  flux  functions  [15],  and  some  nonlinear  subsurface  flow  equations  of 
a  form  different  from  (1.2)  [2,34,36,37,41]. 

Langlo  and  Espedal  [22, 23]  presented  a  macrodispersion  approach  for  the  stochas¬ 
tic  version  of  (1.2).  The  flux  function  is  expanded  in  a  Taylor  series,  and  high- 
order  terms  are  neglected;  then  standard  techniques  represent  macrodispersivity  as  a 
function  of  flow  velocity  covariance.  Zhang,  Tchelepi,  and  Li  take  advantage  of  the 
steady  velocity  field,  and  transform  2-D  flow  to  1-D  Lagrangian  flow  along  stream¬ 
lines  [39,40].  They  formulate  integral  equations  for  moments  from  ensemble  averages 
over  the  streamlines,  rather  than  a  system  of  MDEs.  Variations  on  the  streamline 
approach  are  currently  popular  in  subsurface  transport  [6,8,9,31,32]. 

An  Eulerian  MDE  approach  has  been  successful  for  single-  and  multiphase  pres¬ 
sure  and  velocity  equations  [36],  and  a  natural  next  step  is  to  extend  the  theory  from 
flux  equations  to  transport  equations.  This  framework  differs  from  streamlines  not 
only  in  formulation,  but  also  in  that  the  MDEs  need  no  velocity-distribution  assump¬ 
tion,  and  an  extension  to  transient  velocity  fields  is  relatively  straightforward.  The 
approach  applies  to  any  probability  distribution  of  geologic  properties  and  any  corre¬ 
lation  function,  and  does  not  require  stationarity.  Other  stochastic  theories  generally 
require  such  restrictions. 

Equations  are  derived  in  §2.  In  §3,  we  reduce  the  mean  and  variance  equations 
to  a  simple  form.  These  are  shown  to  be  strictly  hyperbolic  for  1-D  flow,  and  an 
analytical  solution  is  given  in  §4.  The  uniqueness  of  solutions  to  hyperbolic  PDEs  is 
considered.  Classification  is  briefly  discussed  for  2-D  equations.  We  conclude  with 
an  evaluation  of  our  results  and  their  practical  implications,  and  a  brief  overview  of 
additional  questions  that  warrant  future  investigation. 
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2.  Moment  Equations.  We  compare  two  different  approaches  for  statistical 
MDEs  of  1-D  flow.  In  §2.1,  we  expand  fields  in  the  form  u  =  (u)  +  5u  and  the 
resulting  MDEs  are  closed  by  neglecting  products  of  ‘S’  terms.  Analogous  equations 
in  §2.3  result  from  a  full  asymptotic  expansion.  We  show  in  §4  that  the  latter  yields 
moments  that  violate  physical  constraints.  In  both  cases,  random  fluctuations  in 
Y  =  In  AT  are  assumed  small:  ay  < Cl.  We  can  immediately  generalize  to  higher 
dimensions  via  vector  notation.  Moments  of  h  and  v  are  assumed  known  from  (1.1) 
and  moments  of  Y . 

The  1-D  saturation  equation  (1.2),  with  initial  data  s(x,  0)  =  g(x),  is 

dts +  dx(f(s)v)  =  0,  (2.1) 

We  assume  that  g  is  known  with  certainty.  Solutions  are  defined  in  terms  of  vanishing 
viscosity  as  in  §1.1;  henceforth,  this  is  tacitly  understood.  Incidentally,  (2.1)  can  be 
solved  exactly,  and  moments  computed  by  integrating  against  a  known  probability 
density.  This  is  carried  out  in  [39]  and  [40],  for  example.  However,  since  we  ulti¬ 
mately  develop  Eulerian  equations  for  moments  in  two  and  three  space  dimensions 
with  spatial  dependence  in  the  velocity,  we  apply  the  approach  to  the  1-D  equation 
first. 


2.1.  Two-term  expansion.  Moment  equations  may  be  derived  in  a  number  of 
ways.  For  examples  of  commonly  used  methods  applied  to  various  models  in  sub¬ 
surface  flow  and  transport,  see  [7, 11, 14, 15,30,36].  These  methods  originated  in  the 
correlation  equations  of  turbulence  models.  Here  we  apply  a  standard  approach,  sep¬ 
arating  mean  fields  from  random  fluctuations. 

Let  (•)  denote  the  expectation  operator,  defined  by 

W  =  [  V’HdP(w)  (2.2) 

Jn 

for  any  integrable  function  ip  :  SI  ->  M.  on  the  sample  space  SI  with  probability  measure 
P.  We  omit  reference  to  w  in  what  follows. 

The  random  field  Y  is  decomposed  into  deterministic  mean  plus  random  fluctu¬ 
ation:  Y  =  ( Y )  +  5Y.  Each  field  dependent  on  Y  is  represented  similarly: 

h(x )  =  ( h )  (x)  +  5h(x ),  v(x)  =  (v)  (x)  +  5v(x),  (2-3) 

s(x,  t)  =  (s)  (x,  t)  +  5s(x,  t). 

Recall  that  we  only  need  the  decompositions  of  v  and  s  here.  Next,  the  fractional 
flow  function  is  expanded  in  a  Taylor  series  around  ( s ): 

/(*)  =  /«*» + nms + \fnmsg 2 + •  •  •  (2.4) 

So  far,  we  make  no  assumption  regarding  the  size  of  5s  relative  to  ( s ). 

To  obtain  the  mean-saturation  equation,  apply  the  operator  (•)  to  (2.1): 


=  0. 


dt  (s)  +  dx  [/((«))  (v)  +  /'((«))  {5s  5v)  +  i  f"({s ))  {v)  (5s2) 

+  \f"m  (5s2 5v)  +  !/"'«*»  (v)  (5s3)  +  •  •  •  ] 


(2.5) 
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The  flux  terms  include  a  nonlinear  advective  mean,  two  covariances,  and  higher-order 
moments.  For  the  saturation-fluctuation  equation,  subtract  (2.5)  from  (2.1): 

dt6s  +  dx  [/((s))<Jn  +  }'((s))5s  { v )  +  f'((s))  (5s  5v  -  (5s  5v)) 

+\f  "((*))  (v)  ( 5s 2  -  (5s2))  +  ±f"((s))  (5s2 5v  -  (5s2 5v))  +•••]=  0.  (2.6) 

We  derive  equations  for  the  unknown  covariance  functions  (5s5v)  and  (5s2),  using  the 
following  additional  notation.  The  independent  variables  are  x  and  t  except  where 
noted,  and  -\y  denotes  the  replacement  of  x  by  some  y  different  from  x.  It  is  convenient 
and  useful  to  derive  equations  for  the  more  general  two-point  covariances  (dsdflj,)  and 
(ds^s^)  rather  than  for  one-point  covariances. 

To  obtain  an  equation  for  the  saturation-velocity  covariance  (<5s<5ti|y),  multiply 
(2.6)  by  5v(y)  and  apply  (•).  This  results  in: 

d t  +  d x  [/((«))  (5v  5v\v)  +  /'((«))  (v)  (5s  5v\v)  +  f'((s))  ( 5s5v  5v\y) 

+  \f"((s))(v)(5s25v  |„)  +^/"((s))(<5s2<5n5n|J/)-(----]  =0.  (2.7) 

Similarly,  multiplying  (2.6)  by  5s(y,t),  and  using  the  identity2  dt.(5s5s\y)  = 
+  <5s|y<9j<5s,  yields  this  equation  for  the  two-point  saturation  covariance: 

dt  (**!„)  +  dx  [/((«))  (5s|j,5u)  +  /'((«))  (v)  (5s  5s|y)  +  /'((«))  (5s  5v  5s|j,) 

+  \n<s))  (v)  (Ss25s\v)  +  if"((s))  (5s25v5s\v)  +  •  ■  •]  (2.8) 

+5j/[/««Ij/»  (Ss  H  y)  +  f'((s\v ))  Hj/)  +  f'((s\y))  ((5s  <5v)|y<5s> 

+\f"((»\v))  H»)  +  \f"((s\y))  ((5s25v)\y5s)  +  •••]=  0. 

2.2.  Closure  by  perturbation  argument.  If  ay  <C  1  so  that  fluctuations  and 
their  derivatives  may  be  assumed  small  relative  to  the  means,  and  if  /  is  smooth,  then 
we  can  approximate  (2.5)-(2.8)  by  a  closed,  coupled  system.  Defining  csv(x,y,t )  = 
(dsdnl^) ,  cs  =  (fo&sly) ,  cv  =  (5vdv\y) ,  (s)  \y  =  (s)  (y,  t ),  and  csv(x,  y,  t)  =  csv(y,  x,  t), 
the  resulting  system  is 

dt  (s)  +  dx  [/((«))  (v)  +  f'((s))asv  +  i  f"((s))a 2  (v)  ]  =  0,  (2.9a) 

dt.Csv  +  dx  [/({s))c„  +  /'((«))  (v)  csv]  =  0,  (2.9b) 

dtcs  +  dx  [f((s))csv  +  f((s))  (v)  cs]  +  dv  [f((s)\y)c8V  +  f'((s)\y)  (n^)  cs]  =  0.  (2.9c) 

Initial  data  are  (s)(a:,0)  =  g(x ),  csv(x,y,  0)  =  cs(x,y,  0)  =  0;  recall  that  (v)  and 
(5v5v\y)  are  assumed  known.  Both  (2.9b)  and  (2.9c)  have  advective  flux  terms,  are 
coupled  to  the  mean  equation  (2.9a),  and  are  first-order  in  a\.  This  is  consistent  with 
the  approximation  to  (2.9a),  which  is  second-order  in  ay. 

2The  identity  is  not  valid  in  a  strong  sense  for  discontinuous  solutions.  Recall,  however,  that  we 
define  solutions  in  terms  of  the  (smooth)  viscous  solution,  in  the  limit  en  — >  0. 
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Another  common  closure  argument,  which  may  be  called  a  Gaussian  assumption, 
might  be  applied  here.  For  example,  for  the  linear  case  (f(s)  =  s),  (2.9)  is  exact  if  one 
assumes  that  velocity  and  saturation  are  jointly  multivariate  normal  [11].  This  follows 
without  a  perturbation  argument.  The  linear  case  under  this  assumption  is  studied 
in  [15],  with  the  inclusion  of  small-scale  diffusion.  Using  simulations,  we  found  that 
the  Gaussian  assumption,  even  for  transformations  of  saturation  and  velocity  fields, 
is  inappropriate  for  this  problem  [18].  That  result  also  prohibits  the  convenience  of 
modeling  a  transformation  of  f(s)  as  Gaussian  rather  than  using  a  Taylor  expansion 
(a  similar  approach  was  used  successfully  for  unsaturated  flow  in  [2]). 

Note  that  the  mean  equation  (2.9a)  contains  the  functions  asv  =  csv(x,x,t)  and 
a2  =  cs(x,x,t)  rather  than  csv(x,y,t)  and  cs(x,y,t).  This  mix  of  one-point  and 
two-point  covariance  functions  prevents  us  from  immediately  treating  the  MDEs  as 
classically  hyperbolic,  even  though  they  can  be  put  in  conservation-law  form.  Also, 
independent  variables  x  and  y  are  permuted  in  (s)  and  csv  in  (2.9c).  In  general, 
(s)  | y  (s),  and  csv  ^  csv.  We  address  these  issues  in  §3. 

We  refer  to  this  problem  as  “1-D,”  even  though  the  covariance  functions  involve 
two  spatial  variables,  and  the  saturation  covariance  has  fluxes  in  both  directions.  The 
two  variables  represent  two  different  points  in  the  same  1-D  domain.  Similarly,  the 
“2-D”  problem  has  four  spatial  coordinates. 

To  define  the  vanishing-viscosity  solution,  first  formally  add  the  term  endxs  to 
(2.1)  and  carry  out  the  expansion  and  derivation  above.  This  adds  diffusion  terms 
e£><92s,  tDdxcsv,  eDd1cs  to  the  right-hand  sides  of  equations  (2.9a),  (2.9b),  and  (2.9c), 
respectively.  Then  find  solutions  in  the  limit  en  -A  0.  In  practice  we  solve  the  system 
(2.9)  directly. 

2.3.  Infinite  expansion.  In  this  alternative  derivation,  flow  variables  head 
h(x),  velocity  v(x)  and  saturation  s(x,t )  are  represented  by  formal  infinite  perturba¬ 
tion  series  expansions  in  powers  of  a  parameter  e: 

oo  oo  oo 

h  =  ^2enhn(x),  v  =  ^envn(x),  s  =  ^  ensn(x,  t).  (2.10) 

n= 0  n= 0  n= 0 

The  expansion  parameter  e  =  ay  is  shown  to  be  appropriate  within  the  context  of 
the  velocity  and  head  equations  (1.1)  [7,  pp.  184-190],  [38].  For  example,  for  single 
phase,  stationary  uniform  mean  flow  in  1-D,  a 2  is  approximated  by  e2  (y2)  =  Vgtjy. 
Moment  equations  analogous  to  (2.9)  are  derived  in  §A  of  the  Appendix: 

dts0  +  dx  [/(so)uo]  =  0,  dt  (si)  +  dx  [u0/'(s0)  (si)  ]  =  0,  (2.11a) 

d t  <s2)  +  dx  [/(s0)  fa)  +  f'{s0)asv  +  /'(s0)  (s2)  u0  +  ^/"(s0)a2u0]  =  0,  (2.11b) 

dtCsu  +  &x  [/(so)c„  +  /  (so^oCst,]  —  0,  (2.11c) 

dtcs  +  dx  [/(so)cst,  +  /'(s0)wcs]  +  dy  [/(  50 1  y)Csv  +  /,(so||>o|j/C«]  =  0.  (2. lid) 

At  t  =  0,  s o(x,  0)  =  g(x),  and  (si)  (x,  0),  (s2)  (x,  0),  csv(x,  y,  0)  and  cs(x,  y,  0)  are  all 

zero.  The  second-order  mean  is  so  +  e  (si)  +e2  (s2).  Note  that  the  argument  of  /  and 
its  derivatives  is  s0>  the  zeroth-order  mean.  The  system  is  in  fact  closed  again  using 
a  perturbation  argument,  but  now  this  argument  is  contained  in  the  assumption  that 
the  formal  power  series  in  e  converges.  Thus,  throughout  §2,  second-order  equations 
are  closed  by  assuming  that  heterogeneity  is  weak  (ay  <Cl). 
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3.  Classification  of  reduced  equations.  Classification  is  of  interest  for  several 
reasons.  Notably,  macrodispersion  theories  produce  a  parabolic  equation  for  the  mean. 
In  contrast,  our  approach  produces  hyperbolic  equations  for  mean  and  variance.  To 
achieve  this  result,  we  exploit  special  structure  in  the  equations  in  1-D  to  simplify 
(2.9)  and  (2.11).  First,  we  need  the  following 

Definition  3.1.  Letu(x,y,t)  :  ffi2  xffi_|_  -a  ffi",  F  :  ffi”  ->•  ffi",  and  G  :  ffi”  ->  ffi". 
A  system  of  equations  in  conservation-law  form  is  given  by 


dtu  +  9xF(u)  +  53/G(u)  =  0,  u(x,  0)  =  u0(x).  (3.1) 

This  system  is  hyperbolic  if  the  eigenvalues  of  ci  DF( u)  +  C2  DG( u)  are  real  for  all 
ci,  C2  S  ffi,  where  DF  and  DG  are  the  Jacobian  matrices  of  F  and  G  [26]. 

We  expect  the  equations  to  be  nearly  hyperbolic  since  the  original  deterministic 
PDE  (1.2)  is  hyperbolic,  but  the  systems  (2.9)  or  (2.11)  as  given  cannot  even  be  writ¬ 
ten  in  conservation-law  form  due  to  the  inconsistencies  mentioned  above.  However,  we 
show  that  both  sets  of  MDEs  reduce  to  hyperbolic  systems  of  conservation  laws  and 
yield  analytic  solutions.  The  key  discovery  that  allows  this  reduction  is  a  relationship 
between  covariance  functions  in  1-D. 

Observe  that  velocity  v  is  constant  in  1-D,  so  it  is  merely  a  random  variable 
rather  than  a  random  field.  Thus,  velocities  at  any  two  points  in  space  are  perfectly 
correlated  (they  are  the  same  random  variable) .  Equivalently,  the  velocity  correlation 
length  is  infinite.  This  infinite  correlation  length  characterizes  the  principal  difference 
between  stochastic  subsurface  flow  in  one  and  two  space  dimensions  [18]. 

3.1.  Two-term  expansion.  The  MDEs  (2.9)  reduce  to  a  system  of  hyperbolic 
PDEs.  The  first  step  to  this  end  addresses  the  inconsistencies  mentioned  earlier. 

Both  expansion  terms  ( v )  and  5v  are  constant:  by  applying  the  expectation  op¬ 
erator  to  dxv  =  0  we  obtain  dx  ( v )  =  0,  so  that 


0  =  dxv  =  dx  ( v )  +  dx5v  =>  dxSv  =  0. 


(3.2) 


This  implies  that  cv(x,y)  is  also  constant,  and  that  csv(x,y,t)  is  independent  of  its 
second  argument.  Consequently,  cv  is  identical  to  the  second-order  approximation  to 
velocity  variance  a2,  and  csv(x,y,t)  is  identical  to  asv(x,t). 

This  last  identity  removes  the  inconsistency  of  having  asv  instead  of  csv.  We 
still  have  a2  in  (2.9a),  instead  of  cs.  and  we  have  ( s )  \y  and  csv  in  (2.9c).  A  key 
variance-covariance  relationship,  asv  =  asav,  follows:  in  1-D,  the  saturation  profile 
is  completely  determined  by  knowledge  of  the  velocity,  for  any  positive  time.  Thus, 
saturation  and  velocity  are  perfectly  correlated.  A  generalization  of  this  result  is 
obtained  directly  from  the  MDEs  in  §3.3. 

We  divide  (2.9b)  by  av  >  0,  and  retain  only  the  first  two  equations  in  (2.9)  in  the 
following.  Replace  cv  by  a2  and  csv  by  avas,  to  reduce  the  system  to  the  following 
new  equations  for  mean  and  standard  deviation  of  saturation: 


dt 


(s) 


,  q  (  iv)  f({s))  +  ^vf'{{s))as  +  ^{v)  f"{{s))a2 
+  X\  vvf({s))  +  {v)f'({s))as 


=  0.  (3.3) 


Dependence  on  the  second  space  variable  y  has  been  eliminated.  Thus,  (3.3)  is  in 
conservation-law  form,  with  u (x,t)  =  (( s )  ,as),  and  flux  function 

F(n)  _  (  { v )  /«*))  +  o\ ,f'((s))<7,  +  |  (v)  f"({s))a2 
V  avf({s))  +  (v)f'({s))as 


(3.4) 
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Lemma  1.  The  MDEs  (2.9)  reduce  to  a  hyperbolic  system. 

Proof.  It  remains  to  show  that  (3.3)  is  strictly  hyperbolic.  The  Jacobian  matrix 

DF((s)  ,os)  =  (  j.11  j}2  ]  has  entries 

\  312  J22  J 

in  =  (v)  f'((s))  +  avf'((s))as  +  1  (v)  f^((s))o2,  (3.5) 

il2  =  CTvf'{{s))  +  (v)  f"({s))as,  j‘2‘2  =  ( V )  /'((«)). 

DF  is  symmetric;  thus  the  eigenvalues  of  DF  are  real  and  (3.3)  is  hyperbolic.  In 
fact,  we  show  that  (3.3)  is  strictly  hyperbolic;  i.e.,  DF  has  a  complete  set  of  linearly 
independent  eigenvectors  [26].  Distinct  eigenvalues  must  correspond  to  independent 
eigenvectors.  Thus,  we  need  only  consider  the  case  where  DF  has  just  one  eigenvalue. 
This  occurs  when  the  discriminant  of  the  characteristic  polynomial  vanishes:  (jn  — 
J22)2  +  4ji2  =  0.  But  this  would  imply  that  j\2  =  0,  hence  DF  is  a  multiple  of  the 
identity,  and  again  has  a  complete  set  of  eigenvectors.  The  conclusion  follows.  □ 
In  addition,  the  eigenvalues  are  distinct  unless  (( s )  ,os)  €  {(0,0),  (1,0)}  [18]. 

3.2.  Infinite  expansion.  The  evolution  of  moments  in  this  case  is  given  by  the 
system  (2.11).  If  velocity  is  constant  then  Vj  is  constant  for  each  j,  so  long  as  the 
perturbation  series  can  be  differentiated  termwise.  Thus,  cv(x,y)  is  also  constant, 
and  csv(x,y,t )  is  independent  of  its  second  argument.  Consequently,  cv  is  identical 
to  the  second-order  approximation  to  velocity  variance  of,  and  csv(x,  y,  t)  is  identical 
to  osv(x,t).  Finally,  the  relationship  osv  =  ovos  holds  as  before.  Keeping  only  the 
first  four  equations  in  (2.11),  replacing  cv  by  a2  and  csv  by  ovos,  we  obtain  these 
equations  analogous  to  (3.3): 


/  so  \ 

e  (si) 

e2  <S2> 

\  Os  J 


+  dx 


v0f(so) 

v0f'(s0)e(si) 

M  /(«o)  +  <?v  f'(s0)os  +  v0  f'(s 0)  e2  <s2) 


V 


+  \vo  f"(so)c 


Tvf(s0)  +u0  f'(s0)c 


=  0.  (3.6) 


J 


We  show  that  this  system  is  (not  strictly)  hyperbolic.  The  Jacobian  matrix 

(  d  0  0  0  \ 

J21  d  00 

J31  0  d  j'34 

V  J34  0  0  d  j 


DF  (s0,  e(si)  ,  e2  (s2) ,  <rs)  = 


has  entries 


d  =  v0  f'(so),  J21  =  v0  /"(so)e(si) ,  J34  =  ^/'(s0)  +  v0  f"(s0)o3,  (3.7) 

J31  =  e2  {V2)  /'(s0 )  +  avf”(so)os  +  v0  /"(s0)e2  (s2)  +  ^0  f{S\s0)o2s. 

The  spectrum  of  DF  consists  of  the  single  real  eigenvalue  d.  Thus,  the  system  (3.6) 
is  hyperbolic,  but  not  strictly  hyperbolic.  Furthermore,  DF  in  general  does  not  have 
a  full  set  of  linearly  independent  eigenvectors  (this  occurs  only  when  all  off-diagonal 
elements  are  zero).  This  degeneracy  of  DF  leads  to  secular  terms  in  the  solution, 
which  results  in  non-physical  solutions.  We  show  this  explicitly  in  the  next  section. 
Use  of  the  second-order  mean  as  the  argument  of  /(fc)  in  §2.3  yields  a  modified  infinite 
expansion  that  does  not  have  this  drawback,  and  is  a  fourth-order  correction  to  the 
two-term  expansion  [18].  We  will  present  a  comparison  of  the  modified  expansion  to 
the  infinite  and  two-term  expansions,  and  further  analysis,  in  [19]. 
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3.3.  Uniqueness,  and  an  additional  analytical  result.  Because  (2.9)  and 
(2.11)  are  nearly  hyperbolic  systems,  one  might  expect  to  extend  uniqueness  methods 
from  the  theory  of  such  systems.  The  viscosity  method  is  an  appealing  way  to  prove 
uniqueness  [33].  We  briefly  outline  the  method  below,  which  uses  the  vanishing- 
viscosity  solution  introduced  in  §1.1.  Diffusive  effects  are  generally  present  in  the 
physical  problem  that  is  modeled  by  (3.1).  Sharp  profiles  that  evolve  due  to  nonlinear 
advection  are  smoothed  by  diffusion,  and  solutions  remain  smooth  for  t  >  0.  Let  ue 
be  the  solution  to  the  viscous  form  of  (3.1)  with  diffusion  coefficient  (or  “viscosity”) 
en  >  0.  The  viscous  equation  is 

diU-t-d^u)  +dyG(u)  =  eD{d%u  +  dyU).  (3.8) 

It  is  generally  thought  to  be  easier  to  prove  uniqueness  for  solutions  to  (3.8)  than  to 
(3.1)  due  to  the  regularizing  effect  of  the  diffusion  operator  ££>(9^  +  dy). 

We  seek  the  solution  ue  to  (3.8)  in  the  limit  ep  — >  0.  If  this  limit  exists  and  is 
unique  in  a  suitable  function  space,  we  have  the  necessary  result.  Thus,  we  would  like 
to  carry  out  the  following  program: 

1.  Prove  that  (3.8)  has  a  unique  solution  u£. 

2.  Prove  that  ue  has  a  unique  limit  in  a  suitably  defined  function  space  as 
tD  — t  0.  Define  the  solution  to  (3.1)  to  be  this  limit. 

Applying  this  or  similar  approaches  to  the  problem  of  uniqueness  for  conservation 
law  equations  is  found  to  be  exceedingly  difficult.  It  is  usually  applied  only  to  Riemann 
initial  data: 

ju H  X  <  0, 
u  =  < 

[ur,  X  >  0. 

The  states  u/  and  ur  are  usually  assumed  to  be  “close”  in  some  norm.  The  Riemann 
problem  that  consists  of  a  hyperbolic  PDE  in  conservation-law  form  with  Riemann 
initial  data  is  the  fundamental  Cauchy  problem  for  this  class  of  equations  [33]. 

Variations  on  the  viscosity  method  have  been  used  for  1-D  conservation  laws  with 
some  success.  The  result  for  scalar  equations  was  obtained  by  Oleinik  (cited  in  [33]), 
and  was  extended  to  some  systems  of  two  equations  (see  [5]  for  a  recent  result  and  ad¬ 
ditional  references).  Uniqueness  proofs  for  systems  of  equations  are  generally  limited 
to  genuinely  nonlinear  or  linearly  degenerate  equations  [33] .  Neither  the  deterministic 
version  of  equation  (2.1)  nor  the  reduced  systems  above  possess  either  of  these  prop¬ 
erties.  A  review  of  the  literature  does  not  reveal  a  result  general  enough  to  guarantee 
uniqueness  for  (2.9).  Thus  uniqueness  remains  an  open  question;  however,  physical 
and  mathematical  arguments  suggest  that  such  results  can  eventually  be  obtained. 
For  now,  we  must  be  satisfied  with 

Conjecture  3.1  (uniqueness).  Moment  equations  (2.9)  or  (2.11)  have  at  most 
one  vanishing-viscosity  solution  for  uniformly  bounded,  measurable  initial  data. 

Now  we  obtain  a  more  general  form  of  the  covariance  relationship  stated  in  §§3.1 
3.2  directly  from  the  MDEs. 

Lemma  2.  If  there  exists  a  unique  solution  to  (2.9)  with  bounded,  measurable 
initial  data,  then 

cs  cv  —  csv  csv.  (3.9) 


The  same  holds  for  (2.11). 

Proof.  We  prove  (3.9)  for  (2.9)  in  the  sense  of  vanishing  viscosity.  The  proof  for 
(2.11)  is  nearly  identical.  Adding  a  linear  diffusion  term  —dx(endxs)  to  the  left-hand 
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side  of  (2.1)  with  tu  >  0  deterministic,  we  obtain  the  viscous  form  of  the  covariance 
equations 

&tCsv  ~b  ~b  /  ((^))^0 Csv  —  0, 

Cs  ~b  ^/((^))^st>  "b  /  {{s))vqC8 
“b  dy  ^/((s)  \y)csv  +  /  ((s)  ly^olj/Cs  —  —  0. 

Once  the  identity  is  proven  for  (3.10)  and  (3.11),  we  formally  let  tu  — >  0. 

Multiply  (3.10)  by  csv  and  add  to  csv  times  the  corresponding  equation  for  csv. 
Apply  the  identity  csv  dt.csv  +  csv  dtcsv  =  dt(csv  csv).  Note  that  csv  is  independent  of 
x.  We  obtain 

dt{cSv  Csv)  "b  0X  ^/((^))^  csv  -b  /  ((s))vo(csv  csv)  cjjdx{csv  cs^)J 

~b ^/((s)  |j/)c„  csv  +  /  ((s)  \y)v^\yicsv  csv)  —  £d9x{csv  cs^)j  —  0.  (3.12) 

This  equation  is  identical  to  (3.11),  scaled  by  cv.  Both  cvcs  and  csvcsv  are  initially 
zero.  The  identity  (3.9)  follows  by  letting  ep  — >  0  and  recalling  the  uniqueness 
hypothesis.  □ 

Again  asv  =  as  av  follows  by  letting  y  x  so  that  cs  —f  a2s.  and  csv  — >  asv .  It 
is  important  to  recognize  that  this  covariance  relationship  follows  from  observation. 
The  fact  that  it  is  also  a  consequence  of  the  MDEs  (subject  to  uniqueness)  shows  that 
the  MDEs  are  consistent  with  this  intuitive  result. 


(3.10) 

(3.11) 


4.  Solution.  We  present  results  for  the  linear  advection  case  first,  then  a  solution 
for  nonlinear  advection  with  the  two-term  expansion. 

4.1.  Linear  advection.  We  solve  (3.3)  and  (3.6)  for  /(s)  =  s  with  determin¬ 
istic  initial  data  s(x,  0)  =  g(x )  €  C1.  These  easily  obtained  results  exhibit  bimodal 
behavior  similar  to  the  nonlinear  case.  This  linear  flux  case  represents  pure  advection 
of  conservative  solute  transport  in  single-phase  flow,  with  s  as  solute  concentration. 
We  also  show  that  moments  given  by  the  infinite  expansion  violate  physical  bounds. 

4.1.1.  Two-term  expansion.  With  f(s)  =  s,  (3.3)  becomes 


dt 


(s) 

0's 


(v) 

O  V 


(4.1) 


with  initial  data  (s)  (x,0)  =  g{x)  and  crs(x,  0)  =  0.  The  eigenvalues  of  the  Jacobian 
matrix,  representing  wave  speeds,  are  Ai  =  ( v )  —  av  and  A2  =  (v)  +  av.  The  solution 
is  given  by 

{ s){x,t )  =  ^  [g(x  -X2t)+g(x-Xlt)] ,  as(x,t)  =  ^  [g(x-X2t)-g(x-X1t)] .  (4.2) 

Each  moment  is  a  superposition  of  waves  that  move  at  distinct  speeds  Ai  <  A2.  Figure 
4.1  shows  mean  saturation  for  a  Gaussian  initial  profile.  The  peaks  of  the  two  modes 
separate  at  a  speed  of  A2  —  Ai  =  2 av  t.  We  used  the  following  standard  results  for 
uniform  mean  flow,  x  G  [0,  L\.  with  stationary  log  hydraulic  conductivity  (( Y )  and 
ay  do  not  depend  on  x)  [7, 14]: 


(v(x)) 


2 

aY, 


2 


(4.3) 
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where  KG  =  exp((Y))  and  J  is  the  negative  mean  head  gradient.  We  used  parameters 
KGJ  =  0.5,  <f>  =  0.2,  and  aY  =  0.5. 

This  is  a  non-physical  result;  we  do  not  expect  an  initial  localized  solute  pulse 
to  evolve  into  a  bimodal  mean  profile.  Adding  a  linear  diffusion  term  to  (4.1)  will 
not  eliminate  this  bimodality;  in  this  situation,  separation  between  modes  will  not 
be  observed  immediately,  but  will  eventually  appear  due  to  a  separation  speed  that 
is  linear  in  time.  This  is  evident  in  the  nonlinear  case  (§4.3),  where  numerical  and 
artificial  diffusion  are  present  in  our  numerical  scheme,  but  the  two- wave  character  is 
clearly  seen. 

4.1.2.  Infinite  expansion.  Moment  equations  (3.6)  with  f(s )  =  s  are 


f  S  o  \  f  V0  o  0  \  (  s  o  \ 

dt  e2  <«2>  +  e2  (V2)  vo  a,  J  dx  I  e2  <s2)  =  0.  (4.4) 

\  crs  J  \  <JV  0  v0  J  \  as  J 

The  equation  for  (si)  is  not  needed,  since  (si)  =  0.  Each  equation  can  be  solved 
separately,  in  the  sequence  So;  17 s-,  e2  ( s2 )  to  obtain 

s o(x,  t)  =  g(x  -  v0t),  as(x,  t)  =  -avtg'(x  -  v0t),  (4.5) 

e2  (s2)  (x,  t)  =  -e2  (v2)  tg'(x  -  v0t)  +  ^-t2g"(x  -  v0 1). 


The  prime  (')  indicates  a  derivative  with  respect  to  the  argument  z  =  x  —  vq t  of 
the  function  g(z).  The  second-order  approximation  to  mean  saturation  is  given  by 

(s)  «  s  =  s0  +  e2  (s 2). 

2 

The  term  ^-t2 g" (x  —  v0t)  will  dominate  this  approximation  to  (s)  for  large  time,  if 
g"  is  not  identically  zero.  Thus  the  magnitude  of  the  approximation  can  grow  without 
bound.  An  example  for  a  Gaussian  initial  profile  is  shown  in  figure  4.1.  Here  again, 
we  use  a  standard  result  for  uniform  mean  flow  with  stationary  log  conductivity  [38] 
giving  Vo  =  Kg  J/<j>,  e2  (v2)  =  —Vo(Ty/2,  and  av  =  Vq  ay,  and  parameters  are  the 
same  as  in  §4.1.1.  The  mean  behavior  is  similar  for  two-term  and  infinite  expansions; 
the  result  in  the  latter  is  evidently  worse  because  it  violates  physical  bounds  on 
mean  saturation  (saturation  must  remain  between  0  and  1).  This  renders  the  infinite 
expansion  inappropriate.  It  is  clear  from  (4.2)  that  the  mean  always  lies  within 
physical  bounds  for  the  two-term  expansion.  For  the  nonlinear  flux  case,  we  study 
the  two-term  expansion  only. 

Remark  4.1.  Bimodal  mean  concentration  (or  saturation,  in  our  case)  was 
noted  in  [21],  [24]  and  [25].  All  used  methods  to  derive  mean  transport  equations 
that  do  not  involve  second-order  corrections.  In  [21],  Koch  and  Brady  model  the 
concentration-velocity  covariance  using  a  spectral  method,  and  show  that  the  equa¬ 
tion  for  mean  concentration  has  a  wave-like  character  with  two  wave  speeds.  They 
point  out  that  the  bimodality  is  most  pronounced  as  the  velocity  correlation  length 
tends  to  00,  to  which  they  refer  as  the  “most  anomalous  case”.  This  is  the  case  in 
one  dimension,  as  noted  above. 

In  [24,  25],  Lenormand  and  Wang  derive  a  mean  transport  equation  representing 
an  ensemble  of  transport  in  streamlines.  One  version  of  their  equation  is  derived  from 
a  series  expansion  of  a  Fourier  transform,  truncated  at  second  temporal  moments.  For 
an  example  of  layered  material,  this  equation  is  shown  to  have  the  wave  character  as 
in  [21],  and  it  is  suggested  that  bimodality  may  be  eliminated  by  incorporating  higher- 
order  temporal  moments.  If  this  is  indeed  the  case,  it  may  suggest  that  incorporating 
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Fig.  4.1.  Mean  saturation,  linear  flux  with  Gaussian  initial  profile  for  uniform  mean  flow. 
Both  two-term  (solid)  and  infinite  expansion  (dotted)  moments  exhibit  bimodal  behavior.  The  latter 
violates  physical  bounds  on  saturation. 


higher-order  statistical  moments  in  our  equations  eliminates  bimodality.  However, 
when  higher  moments  are  retained  in  the  mean  equation,  equations  for  these  moments 
must  also  be  generated.  The  system  of  equations  expands  very  quickly  as  higher 
moments  are  incorporated. 

4.2.  Nonlinear  advection.  Equation  (3.3)  is  in  the  form  (3.1)  with  G  =  0. 
For  nonlinear  F(u),  a  solution  consists  of  a  sequence  of  shock  waves,  constant  states, 
and  rarefaction  waves.  We  begin  by  briefly  reviewing  some  of  the  theory  of  hyperbolic 
systems  of  conservation  laws.  Several  details  are  omitted  and  can  be  found  in  [29, 33]. 

4.2.1.  Rarefactions.  Rarefaction  waves  are  self-similar  solutions  in  C1  that  de¬ 
pend  on  the  similarity  variable  r]  =  x/t  only.  Let  u  =  u(rj);  then  (3.1)  implies  that 

dvt  „  „ .  ,  du.  , .  , 

-,-+DF(u)-=0,  (4.6) 

where  DF( u)  is  the  Jacobian  matrix  of  F(u).  Let  the  eigenvalues  of  DF( u)  be  ordered 
Ai(u)  <  A2 (u)  and  let  r*,(u)  be  the  associated  eigenvectors  (k  =  1,2).  Then  must 
be  an  eigenvalue  Afc(u)  and  du/drj  is  proportional  to  the  eigenvector  r*,(u).  This 
relationship  is  given  explicitly  as 

(VAjfc(u)-rt(u))  1rk( u)  (4.7) 

for  each  k,  and  is  valid  as  long  as  the  directional  derivative  VAfc(u)  •  r^(u)  does  not 
vanish.  It  follows  that  these  waves  are  integral  curves  of  the  vector  fields  defined  by 
eigenvectors  of  DF  in  the  phase  plane  defined  by  u.  The  integral  curves  must  have 
eigenvalues  increasing  in  iy,  otherwise,  characteristics  cross  and  new  shocks  form. 
The  most  general  theory  for  hyperbolic  systems  requires  that  the  system  be  either 
genuinely  nonlinear  (VA k  ■  rk  ^  0)  or  linearly  degenerate  (VXk  •  =  0)  [4,  29, 33]. 

4.2.2.  Weak  solutions,  shocks,  and  entropy  conditions.  Shocks  are  al¬ 
lowed  if  we  consider  the  weak  form  of  (3.1).  Let  the  shock  be  defined  in  the  (x,  t) 
plane  by  x  =  £(t);  its  speed  is  given  by  £  =  d£/dt.  Let  u_  and  u+  be  the  solution  val¬ 
ues  to  the  left  and  right  (along  the  x-axis)  of  the  shock.  Then  the  Rankine-Hugoniot 
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condition  is  given  by 

[u+-u-]^)  =  [F(u+)-F(u-)].  (4.8) 


Entropy  conditions  are  then  imposed  to  capture  the  vanishing-viscosity  solution 
to  (3.1).  For  scalar  hyperbolic  equations,  this  requirement  is  sufficient:  the  vanishing- 
viscosity  solution  is  identical  to  the  weak  solution  that  satisfies  an  entropy  condition. 
Since  this  issue  is  directly  related  to  the  question  of  uniqueness,  such  a  statement 
cannot  yet  be  made  for  (3.1)  in  general.  We  apply  the  following  well-known  entropy 
condition. 

Definition  4.1  (Lax  condition).  A  shock  connecting  the  states  u_  and  u+, 
traveling  with  speed  £  is  admissible  if 

Afc(u+)<^<Afc(u-)  (4.9) 


for  one  of  the  characteristic  speeds  X/.,  k  =  1,  2. 

In  words,  the  characteristic  curves  defined  by  dx/dt  =  Afc(u)  must  run  into  the 
shock,  not  out  of  it,  as  t  increases  [33].  A  solution  in  the  phase  plane  is  a  connected 
set  of  rarefaction  and  shock  curves  that  satisfy  the  Lax  condition.  However,  in  our 
case  it  is  not  enough  to  simply  apply  (4.6)-(4.9).  As  we  present  our  solution,  we 
introduce  an  additional  entropy  condition  due  to  Liu  [27-29]  to  address  the  specific 
structure  of  our  equations. 

4.3.  Moment  Equations.  Now  we  present  a  solution  to  the  MDEs  (3.3)  for 
nonlinear  f(s)  that  satisfies  entropy  conditions.  For  illustration  we  choose  f(s)  = 
s2 /(s2  +  m(  1  —  s)2)  with  viscosity  ratio  m  =  1/2.  We  consider  deterministic  initial 
saturation  given  by  the  Heaviside  function  s(x,  0)  =  H[—x],  so  that  as  is  initially 
zero.  This  models  an  oil-saturated  reservoir  with  water  forced  in  at  the  origin  (we 
assume  v  >  0).  To  simplify  notation  of  (3.3),  set  u±  =  (s),  U2  =  crs,  and  r  =  (v)  t  so 
that  dt  =  (v)  dT.  Let  e  =  avf  (v)  (this  is  consistent  with  the  e  previously  used  as  an 
expansion  parameter).  The  scaled  equations  are 


Ml 

«2 


+  9X 


f  +  ef'u2+  \f" 
ef  +  fu2 


u\ 


=  0. 


(4.10) 


Note  that  the  argument  of  /  and  its  derivatives  is  always  u\.  Let  u  =  (m,  1*2).  The 
Jacobian  matrix  is 


DF{  u) 


/'  +  e/"u2  +  |/(3)u|  ef  +  f"u2  \ 

ef  +  fnU2  f  J  • 


(4.11) 


^From  §3.1  we  know  that  the  eigenvalues  of  DF  are  real  and  may  be  ordered  Ai  <  A2 
except  at  u  e  {(0,0),  (1,0)},  where  both  are  zero.  The  corresponding  eigenvectors 
are  r*(u),  k  =  1,  2. 

The  true  mean  and  standard  deviation  of  saturation  must  be  non-negative.  The 
mean  cannot  exceed  unity,  as  it  represents  a  fraction.  Therefore,  we  only  consider  the 
subdomain  u  e  [0, 1]  x  [0, 00)  of  the  phase  plane.  The  endpoints  of  the  solution  to 
(4.10)  are  u  =  (1,0)  and  u  =  (0,0),  which  represent  a  boundary  condition,  and  the 
initial  condition  for  positive  x,  respectively. 

Vector  fields  Tk{u(rj))  are  shown  in  figure  4.2,  pointing  in  the  direction  of  increas¬ 
ing  eigenvalue  (VA&  -r*  >  0).  To  obtain  a  rarefaction  curve,  one  solves  (4.7)  for  given 
initial  data  for  a  single  k  value.  Moving  in  the  direction  of  increasing  -q ,  a  rarefaction 


MOMENT  EQUATIONS  FOR  IMMISCIBLE  FLOW 


13 


continues  so  long  as  VA*,  •  >  0.  An  inflection  locus  [17],  consisting  of  points  where 

VAfc  •  rfc  =0,  can  be  inferred  from  figure  4.2  where  vectors  reverse  direction.  This  lo¬ 
cus  of  curves  represents  a  barrier  to  rarefactions.  The  existence  of  the  inflection  locus 
is  evidence  that  the  system  is  neither  genuinely  nonlinear  nor  linearly  degenerate. 


0.8 


'  0.6 


0.4 


0.2 


r  (slow  rarefactions) 


0.5 

Mean  saturation  (u^ 


r2  (fast  rarefactions) 


(  \  \  \  \  _ 

1 

s  /  /  /  i  i  (  j  \  \  \  S^vvvvvvv'' 

_ /  /  /  /  1 1 1 1 1 

imww _ _ _ //////// 

(\\uwv — rrtfttt  t 

0.8 

s  /  /  /  /  l  (  l  f  \  \  N  ^ 

//  1  l  l  \  \  \  \  Ss'vvvv''vvv 

{  1  1  |  I  \  \  \  \ 

W\  \  \\Wv^ _ Sfffftt  t 

lUWWv^ _ _ / /  f  t  /  /  f  f  t 

IWWNVVSV _ _ _ _ _ _  /  /  f  t  f  t 

0.6 

/  /  /  M  j  f  \  \ 

""">//  f  l  I}  \  \  \  \VVVVVVVV'“'‘ 

'2'  A  /  /  f  l  j  1  ^  \  \  ^Wvvvvvvvv 

iWUWw _ ////// / 

^ ' 's  /  /  J  l  l  J  ^  \  \  \  ^'''vv''vvvvv 

lUWUw  / 

(\\\ 

0.4 

'"""'s/  /  /  /  1  '  V  \  ^NNNVVVVVVVV 

/  1  im\\\\\sssvv 

WU'yUW'.N'*-'"*'////////// 

0.2 

"///// /  /  /  / 1 

. WVNWWVWWW 

0 

,  ///////// . 

0.5 

Mean  saturation  (u 


Fig.  4.2.  Slow  and  fast  rarefaction  vector  fields. 


The  steps  to  a  solution  are  similar  to  those  in  the  deterministic  case,  which  is 
reviewed  for  comparison  in  §B.  We  first  study  the  behavior  near  the  boundary  value 
u  =  (1,0),  and  follow  the  solution  forward  in  -q  =  x/t.  When  e  =  0,  this  is  just 
the  scalar  deterministic  case,  which  has  a  rarefaction  connecting  to  s  =  1.  Thus, 
we  expect  a  rarefaction  to  connect  to  (1,0).  We  find  that  only  a  slow  rarefaction, 
corresponding  to  Ai  and  designated  Ri(rj),  may  connect  to  (1,0)  (see  figure  4.3). 
This  curve  is  continued  from  (1,0)  to  a  point  u*  just  short  of  the  inflection  locus. 
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Fig.  4.3.  Construction  of  the  slow  rarefaction  R\  from  (1,0). 


Next,  consider  the  solution  near  (0,0),  the  “end”  value  of  u  on  the  right  (in  x ). 
We  find  that  a  shock  must  connect  to  this  point,  as  in  the  deterministic  case.  The 
point  (0, 0)  lies  to  the  right  (in  x)  of  the  shock,  so  let  u+  =  (u+,  «+)  =  (0,0).  All 
points  u  =  (iii,  ii2 )  on  a  possible  shock  curve  connecting  to  u+  are  defined  by  (4.8) 
as  follows: 
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(u+-«1)4  =  F1(u+)-F1(u)  (4.12) 

(u+  -u2)i  =  F2( u+)  -  F2( u), 


where 

Fi(u)  =  f(u 0  +  e/'(wi)«2  +  ^  (4.13) 

F2(u)  =  e/(ui)  +  f'(ui)u2 

(from  (4.10)).  Eliminating  £  in  (4.12)  gives  an  implicit  definition  of  the  curve  on  which 
points  u  lie,  and  we  may  write  this  definition  formally  as  $(u;  u+)  =  0.  This  gives 
a  continuous  set  of  possible  u,  speeds  £,  and  eigenvalues  A*(u).  We  follow  this  curve 
from  the  starting  point  u+  as  long  as  the  Lax  entropy  condition  (4.9)  is  satisfied, 
which  is  tested  by  computing  the  speed  £  using  (4.12),  and  the  eigenvalues  A*,  at  each 
point  u.  The  Lax  condition  in  this  case  is 

Afe(u+)  <  £  <  Afc(u).  (4.14) 

The  point  u_  to  the  left  (in  x)  of  the  shock  is  the  endpoint  of  this  curve  in  the  phase 
plane.  In  addition,  the  following  criterion  must  be  met  [27,29]. 

Definition  4.2  (Liu  Criterion).  Define  the  Hugoniot  curves 

Hk(u0)  =  {u  :  (u  -  uo)s  =  F(u)  -  F(u0)},  k  =  1,2 

for  some  scalar  ^  =  <;(u0,  u).  The  shock  conditions  (4.14)  tell  us  that  u+  S  Hk(u~) 
and  £  =  ?(u“,u+).  A  shock  with  speed  £  =  <^(u_,  u+)  is  admissible  if  it  satisfies 

^  =  c(u“,u+)  <^(u“,u)  (4.15) 

for  any  u  €  Hk(u~)  between  u_  and  u+. 

In  general,  two  curves  may  pass  through  each  point  u+  and  satisfy  4>(u;  u+)  =  0 
[29, 33],  but  only  one  continues  from  (0,  0)  into  the  first  quadrant.  This  curve  is  shown 
on  the  left  in  figure  4.4  along  with  Ri.  At  each  point  on  the  curve,  Ai  (u)  <  £  <  A2(u), 
so  it  is  a  “2-shock”.  We  designate  this  curve  S2  and  its  speed  £2.  Points  along  this 
curve  are  admissible  as  long  as  (4.14)  and  (4.15)  are  satisfied  for  k  =  2.  Notice  that 
this  curve  crosses  the  inflection  locus  for  A2 . 

Evidently,  an  intermediate  solution  between  Ri  and  S2  is  necessary.  A  slow 
rarefaction  cannot  continue  on  to  meet  S2,  and  we  find  that  a  fast  rarefaction  curve 
cannot  connect  the  two  curves.  The  connection  must  involve  a  shock;  the  simplest 
possibility  is  a  single  slow  shock.  We  follow  a  possible  shock  forward  in  r]  =  x/t  from 
points  on  R\.  This  leads  to  an  admissible  solution,  if  we  account  for  the  following. 
We  know  that  the  intermediate  shock,  designated  S±  with  speed  £l5  must  satisfy  the 
entropy  condition  Ai(u+)  <  <  Ai(u“),  and  the  condition  (4.15).  But  since  Si 

connects  to  a  rarefaction  on  the  left  (in  x),  we  must  have  Ai(u“)  =  £i.  Otherwise, 
the  value  u_  to  the  left  of  the  shock  travels  faster  than  the  shock  front,  which  is  not 
allowed.  Thus,  we  find  the  point  at  which  i?i  connects  to  S±  by  matching  rarefaction 
and  shock  speeds  at  the  front. 

An  algorithm  for  the  speed-matching  step  for  the  MDEs  is  provided  in  §C  of  the 
Appendix.  By  following  a  shooting  procedure  outlined  in  that  algorithm,  we  obtain 
an  approximation  to  a  curve  labelled  Si  which  connects  R\  to  S2 ,  shown  on  the  right 
in  figure  4.4.  The  condition  (4.15)  is  also  satisfied  for  Si  at  the  end  of  this  process. 

The  solution  in  physical  space  is  pieced  together  from  the  phase  space  solution 
as  follows.  Fix  r  =  {v)t.  Along  Ri .  the  location  is  computed  by  x  =  Ai(u)r,  for 
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Fig.  4.4.  Construction  of  shocks  connecting  to  rarefaction  R\  and  (0,0).  First,  a  fast  shock  S2 
connecting  to  (0,0)  is  constructed.  Then  a  slow  shock  is  introduced  to  connect  R\  and  S 2.  “Fast” 
and  “slow”  vector  fields  r2  and  ri  are  illustrated  in  the  left  and  right  plots,  respectively. 


u  between  u_  and  (1,0).  Ri  and  Si  connect  at  u_,  where  x W  e(it  =  Ai(u_)r. 
The  solution  is  constant  from  x^1)  to  x^  =  £2  t,  where  it  then  jumps  to  zero.  This 
solution  is  compared  to  the  solution  obtained  from  our  numerical  PDE  scheme  in 
figure  4.5,  for  L  =  2,  m  =  1/2,  (v)  =  5/2,  av  =  5/4  at  a  fixed  time  t  =  0.2. 

We  have  a  solution  that  satisfies  entropy  conditions,  but  uniqueness  of  the  so¬ 
lution  is  not  resolved.  Most  results,  again,  require  genuine  nonlinearity  or  linear 
degeneracy,  and  our  system  does  not  satisfy  these  conditions.  Liu  extended  existence 
and  uniqueness  results  to  more  general  systems  [27,  28],  but  the  restrictions  he  places 
on  flux  functions  are  not  met  by  Fi  and  F2,  and  most  results  require  “small  data” 
(that  is,  a  small  jump  ||u“  —  u+||;  see  [5]  for  one  extension  to  “large  data”). 

However,  it  is  clear  in  figure  4.5  that  this  solution  matches  that  obtained  from  a 
numerical  PDE  scheme  applied  directly  to  (4.10).  Our  first-order  upwind  PDE  scheme 
includes  numerical  and  artificial  diffusion.  The  diffusion  coefficient  is  roughly  three 
orders  of  magnitude  smaller  than  the  jumps  in  solution  values.  Thus  the  numerical 
solution  is  near  the  vanishing- viscosity  limit.  The  numerical  result  also  shows  that 
the  addition  of  linear  diffusion  terms  does  not  eliminate  bimodality  in  the  solution. 

The  saturation  variance  is  supported  primarily  on  an  uncertainty  interval  between 
fronts.  Physically,  the  solution  represents  two  zones  containing  mixtures  of  the  two 
fluid  phases  (for  example,  water  and  oil),  and  a  third  containing  only  the  oil  phase. 
In  the  first  zone,  we  have  a  smoothly  varying  mixture  from  the  injection  boundary 
(x  =  0),  where  the  mean  oil  content  tends  to  zero,  down  to  a  constant  mixture  just 
left  of  a  shock.  In  the  second  zone,  we  have  a  constant  mix  of  oil  and  water.  The 
solution  does  not  represent  physical  reality.  Rather  than  two  shock  waves,  the  true 
mean  saturation  is  more  likely  to  have  a  smooth,  diffuse  profile.  Taken  with  the  profile 
of  as ,  however,  these  second-order  solutions  provide  some  insight  into  the  propagation 
of  uncertainty  in  two-phase  flow. 

In  the  limit  ay  -A  0  the  mean  saturation  tends  to  the  classic  Buckley-Leverett 
profile  in  figure  B.l.  The  analogous  construction  of  a  solution  in  the  scalar  deter¬ 
ministic  saturation  case  (§B)  can  now  be  stated  succinctly.  In  that  case,  the  initial 
rarefaction  is  followed  forward  in  x,  up  to  the  inflection  point  f"(s*)  =  0.  A  shock 
connects  to  the  rarefaction  from  s+=0to  s_>s*,  where  the  shock  speed  matches 
the  characteristic  speed  v  f(s~)/s~  =  v 
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Fig.  4.5.  Mean  and  standard  deviation  of  saturation.  The  solid  curve  is  obtained  from  semi- 
analytical  construction  of  shocks  and  rarefaction  curves  in  phase  space;  dashed  curve  is  obtained 
from  numerical  PDE  scheme  applied  directly  to  (4.10). 


For  our  semi-analytical  solution  described  above,  the  approximation  to  R\  was 
obtained  with  an  explicit  Runge-Kutta  (4,5)  formula  (function  ode45)  in  MATLAB®. 
Shock  curves  were  approximated  by  incrementing  u±,  and  solving  for  un  using  the 
roots  function  in  MATLAB. 

5.  Conclusions.  For  spatial  resolution  of  uncertainty,  these  bimodal  results  sug¬ 
gest  that  second-order  Eulerian  MDEs  may  be  inappropriate  for  1-D  immiscible  flow. 
Unlike  1-D,  v  has  finite  correlation  length  in  2-D,  and  macrodispersion  models  are 
based  on  correlations  in  Sv.  In  results  for  analogous  MDEs  for  passive  2-D  solute 
transport  with  diffusion,  bimodality  was  not  observed  [15];  one  might  expect  macrodis¬ 
persion  to  lead  to  a  similar  result  for  immiscible  flow.  Nevertheless,  in  a  forthcoming 
submission  we  show  that  somewhat  mitigated  bimodality  does  persist  in  2-D,  even 
with  diffusion  terms  [20].  More  positively,  for  spatial  averages  such  as  oil-production 
curves,  good  matches  to  Monte  Carlo  simulations  are  found. 

6.  Acknowledgments.  Ken  Jarman  would  like  to  thank  Joseph  Oliveira  for 
many  helpful  suggestions  and  invaluable  guidance  through  extensive  conversations. 

Appendix  A.  Moment  equations  for  infinite  expansion. 

We  derive  equations  for  moments  of  saturation  and  velocity  using  the  asymptotic 
series  in  (2.10).  Moments  of  In  A,  h,  and  v  are  known.  A  note  on  terminology  is 
appropriate.  The  term  “order”  applies  to  the  power  of  e,  rather  than  the  order  of  the 
statistical  moment.  In  special  cases  we  may  use  the  term  “order”  in  the  latter  sense 
and  clearly  differentiate  from  order  in  e.  Thus,  a  “second  moment”  is  a  covariance, 
but  “second-order  moment”  is  any  moment  that  is  approximated  to  order  e2. 
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Expand  f(s)  in  a  Taylor  series  around  s0(x,t): 

f(s)  =  /(so)  +  /'(«o)(*  _  s0)  +  ~f"(s0)(s  -  s0)2  +  ■  ■■ 

=  /(s  o)  +  /'(so)  (e(si  —  <*1»  +  e2(s2  —  (S2))  +  e3  S3  +  •  •  • )  (A.l) 

+  2^"(s°He(Sl  —  (si))  +  g2(s2  —  ($2))  +  e3S3  +  •  •  • )  H - • 

Applying  this  expansion  to  (2.1)  obtains  the  zeroth-  and  first-order  equations 

dts0  +  dx[f(a0)v0]  =  0,  (A. 2) 

dtsi  +  dx  [/(s0)ui  +  /'(so)siuoj  =  0.  (A.3) 

Observe  that  (s0)  =  s0  since  all  terms  in  the  equation  for  s0  are  deterministic.  From 
the  theory  of  MDEs  of  single  phase  flow  we  find  that  (iq)  =0  [41],  so  that  applying 
the  operator  (•)  to  (A.3)  gives 

dt  <si)  +dx  [vo/'(so)  <si) ]  =  0,  (A.4) 

with  initial  data  (si)  (x,0)  =  0. 

Proposition  3.  Where  (A. 2)  and  (A.4)  admit  solutions  So,  (si)  inC1,  (si)  =0. 
Proof.  Solutions  in  C1  are  self-similar,  with  similarity  variable  T)  =  x/t.  Substi¬ 
tuting  s0(h)  into  (A. 2)  gives 

v0f'(s0(v))  =  V 

Using  this  identity  and  substituting  (si)  =  (si)  (rj)  into  (A.4)  gives 
0  =  dt  (si)  +dx  fo(si)] 

=  »,(«.)§  + a, 

-  - (si)  +  (si)  +  j  (si}  > 

and  the  proposition  follows.  □ 

We  conjecture  that  (si)  =  0  for  discontinuous  solutions  as  well.  Applying  (•)  to 
the  second-order  equation  gives  the  following: 

dt  <s2)  +  dx  [/(s0)  (^2)  +  /'(s0)  (siiq)  +  f'(s0)  (s2)  u0  +  ^/"(s0)  (si2)  w]  =  0.  (A.5) 

Now  we  derive  equations  for  second  moments.  Multiply  the  first-order  equation  (A.3) 
by  vi  (y)  and  apply  (•)  to  obtain  the  equation  for  the  two-point  saturation- velocity 
covariance, 

3i(siDi|j/)+5a![/(so)(wiWi|j/)+/,(so)wo(siWi|j,)]  =0.  (A.6) 

Similarly,  multiply  (A.3)  by  Si(y,t)  and  apply  (•)  to  obtain  the  equation  for  the  two- 
point  saturation  covariance, 


dt  (sisilj, 

+  d, 


)  +  dx  [/(s0)  (silyiq)  +  f'{s0)v0  (sisilj,)  ] 

V  [/(So|y)  (sisilj/)  +  /  (^0  1 2/)'Co  1 2/  (sl|j/sl)  j 


=  0. 


(A.7) 
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The  system  for  this  version  consists  of  the  moment  equations  (A. 2)  and  (A.4)-(A.7). 
Defining  csv(x,y,t )  =  (sit>i|y),  cs  =  (siSi|y),  cv  =  and  csv(x,y,t)  = 

c8V{y,x,t),  the  system  is  given  by  (2.11).  If  the  conjecture  that  (si)  =  0  holds, 
then  it  is  not  necessary  to  include  (A. 4). 

Appendix  B.  Deterministic  Buckley  Leverett.  When  ay  — t  0,  (2.1)  and 
(3.3)  reduce  to  a  deterministic  Buckley-Leverett  model: 

dts  +  dx(f(s)v)  =  0,  (B.l) 

with  initial  data  s(x,  0)  =  H[—x]  as  in  §4.3.  Recall  f(s)  =  s2 / (s2  +  m(l  —  s)2);  figure 
B.2  shows  this  curve  with  viscosity  ratio  m  =  1/2.  The  well-known  vanishing- viscosity 
solution  to  (B.l)  is  given  by 

s(x,t)  =  w(rj)H[f' (s~)vt  —  x]  with  x  >  0,  t  >  0,  (B.2) 

where  r]  =  x/t,  w(r ])  is  the  solution  to  vf'(w(rj))  =  rj,  and  s~  is  defined  below.  The 
(water)  saturation  profile  at  a  positive  time  is  shown  in  figure  B.l.  Any  monotone 
S-shaped  function  for  f(s)  will  lead  to  a  qualitatively  similar  solution. 


Fig.  B.l.  Buckley-Leverett  saturation  profile  at  t  >  0. 


We  briefly  review  the  solution  of  (B.2)  within  the  context  of  §4.2.  The  rarefaction 
part  of  the  solution  is  obtained  from  (4.6),  here  given  by 

ds  .  ..  ds  _  . 

-,-+»/ („(,»-=0,  (B.3) 

so  that  r]  =  v  f'{s(rj)).  This  simply  states  that  s(x,  t )  is  constant  along  characteristics, 
and  leads  to  the  definition  of  w(rj)  in  (B.2).  For  t  >  0,  rarefaction  is  allowed  as  long  as 
/,(s(7?))  is  increasing  in  rj;  otherwise  a  shock  forms.  The  requirement  VA*  •  rj,  ^  0 
for  (4.7)  reduces  to  v  f"(s)  0,  a  convexity  condition  that  is  violated  at  a  point  s* 

where  /"(«*)  =  0.  The  rarefaction  must  connect  to  a  shock  at  some  value  s~  >  s*. 
The  condition  (4.8),  here  given  by 

[s+  ~  s~]i(t)  =  v  [/(s+)  -  (B.4) 

gives  us  information  about  the  shock  £. 

It  is  evident  from  the  shape  of  /'  (figure  B.2)  that  characteristics  to  the  right  of 
the  shock  must  run  into  it,  so  that  s+  must  be  the  initial  value  zero.  Also,  £  must 
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Fig.  B.2.  The  nonconvex  fractional  flow  function  and  its  first  derivative. 


match  the  characteristic  speed  q  =  v  f'(s~)  just  to  its  left.  Combining  these  facts,  we 
obtain  —s~ =  — f(s~ ),  which  can  be  solved  uniquely  for  s~ .  It  is  the  value  of 
s  at  which  the  secant  of  /  that  passes  through  the  origin  is  identical  to  the  tangent 
of  /.  In  our  case,  s~  >  s*.  These  arguments  are  analogous  to  the  application  of  the 
Liu  criterion  in  §4.3. 

Physically,  the  profile  in  figure  B.l  represents  a  smoothly  varying  mix  of  oil  and 
water  from  the  inflow  boundary  with  a  sharp  transition  to  an  oil-saturated  zone. 

Appendix  C.  Algorithm  for  connecting  i?i  to  £2  by  a  shock. 

Define  tol  to  be  the  error  tolerance. 

1.  Choose  u“  =  u*  on  Ri.  Define  q*  by  u*  =  Ri(q*). 

2.  Solve  4>(u;  u“)  =0.  Define  u+  =  u  at  the  point  of  intersection  with  £2. 

3.  Compute  £1  =  [Fi(u+)  —  Pi  (u_ )]/[■«+  —  uf]. 

(i)  If  16-  A!  (11- )  I  <  tol ,  then  stop. 

(ii)  Else  if  <  Ai(u“),  then  choose  q  <  q*. 

(iii)  Else  choose  q  >  q* . 

4.  Repeat  at  step  2. 
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